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SUMMARY 

An implicit ADI numerical method for the calculation of two-dimensional 
unsteady flows with strong convection effects is described. The method is 
based upon the conventional Crank-Ni chol son approach for parabolic equations 
but an upwind-downwind differencing is used for the first order spatial 
derivatives associated with convection. The differencing is carried out in 
the current and previous time plane in such a way that the algorithm is second 
order accurate in both space and time. The difference equations are factored 
into sequential operators, one in each independent spatial variable; the 
solution at each time step may then be computed as a sequence of tridiagonal 
matrix problems. The method may be used in a noniterative manner although 
iteration at each time step is recommended in situations where the effects of 
convection are strong. 


INTRODUCTION 

Unsteady flows with strong convection effects occur in a variety of 
circumstances. Many of these situations are associated with unsteady 
boundary-layer separation phenomena and the strong unsteady vi scous-invi scid 
interactions that are observed to occur between an outer effectively inviscid 
flow and the viscous flow near a solid surface. Specific examples include (i) 
small separation bubbles on the upper surface of turbine blades and airfoils 
which in certain situations erupt into the inviscid flow region, (ii) the 
eruption of boundary-layer flows which is induced by the motion of vortices 
near solid walls (refs. 1 and 2) and (lit) bursting in turbulent boundary 
layers (refs. 3 and 4). One feature of the aforementioned flow situations is 
that, as the interaction initiates, very strong updrafts begin to develop in 
the boundary layers near the wall and the flow field locally is dominated by 
strong convection effects. 

When finite difference schemes are used to compute the evolution of a 
time-dependent flow, methods based on some version of either the Crank- 
Ni chol son algorithm (ref. 1) or a factored-operator ADI (Alternating- 
Direction-Implicit) technique are used. In Crank-Ni chol son methods, the 
difference equations in the current time plane are usual ly. sol ved by 
iteratively sweeping point-by-point through the two-dimensional spatial mesh 
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until convergence is attained at each point. In the factored alternating- 
direction methods (ref. 5), the difference operators are factored and the 
difference equations are solved as an ordered sequence of tridiagonal matrix 
problems in each of the coordinate directions. In general, the ADI methods 
are more efficient than the Iterative Crank-Ni col son method. However, in a 
recent study of the evolution of unsteady flow in a two-dimensional boundary 
layer (ref. 6), the explicit Beam and Warming (ref. 5) type of algorithm was 
found to fail in situations where strong time-dependent convective effects 
were present. The algorithm described in this study is also a factored ADI 
method; however, upwind-downwind differencing is used for first order 
convective derivatives and the method is also implicit. The present method 
produced good results (ref. 6) in a situation where both standard Crank- 
Nicolson method and conventional ADI methods failed. 


Problem Statement 


The two-dimensional unsteady boundary-layer equations in either Eulerian 
or Lagrangian coordinates can be written in the following general form: 
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where u is the tangential velocity in the boundary layer, (£,n) are 
independent spatial coordinates and t is time. In the Lagrangian 
description of the boundary-layer motion, (£,n> are the Cartesian coordinates 
of a particle at the initiation of the motion; if x(£,n,t) denotes the 
streamwise position of a fluid particle at time t, equation (1) is 
supplemented with an equation 


and the functional coefficients in equation (1) (namely R,S,T,P,Q,W) depend on 
x and u as well as (E,n,t). Consequently the system described by 
equations (1) and (2) is nonlinear and equation (2) is simply a convenient 
representation of the momentum equation. Note that the coefficients R and 
T are such that 


R > 0, T > 0. 


(3) 


In the Eulerian description of unsteady two-dimensional boundary-layer flow 
(see, for example, ref. 1), 


R = S = 0, T = 


and equation (1) is supplemented with 
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where <j/ is the two-dimensional stream function; in this case, the functional 
coefficients in equation (1) are functions of u and as well as (E,n,t). 
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It is worthwhile to note that each momentum equation for the full 
two-dimensional Navi er-Stokes equations is also in the general form of 
equation (1). 

In order to compute the evolution of the flow field, it is necessary to 
develop algorithms which advance the solution of equation (1) one time. step at 
a time. Consider a spatial mesh in the £n plane and adopt a convention 
where subscripts ij denote a quantity evaluated at a typical point in the 
mesh, with i and j representing typical locations in the £ and n 
directions respectively. It is assumed that the solution is known at time 
t* (the previous time plane) and the objective is to compute the solution at 
t = t* + At (the current time plane), where At is the time step. All 
quantities evaluated in the previous time plane are assumed known and will 
indicated by an asterisk; values of u ir the current time plane are unknown 
and to be found. 


The Conventional Crank-Ni col son Method 

The grid structure near a typical spatial mesh point is indicated 
schematically in figure 1. In the standard Crank-Ni col son method (ref. 1), 
equation (1) is approximated at a point midway between the current and 
previous time planes along a line connecting the typical point in mesh 
labelled ij. The coefficients in equation (1) are evaluated through a simple 
averaging procedure involving quantities in the current and previous time 
planes; using an overbar to denote quantities evaluated at the midpoint, a 
typical coefficient in equation (1) is evaluated according to 





( 6 ) 


for example. For derivatives in equation (1), an average is also carried out 
between the current and previous time planes and central difference 
approximations are subsequently used in both time planes. Let and 6^ 

be central difference operators in the i and n directions respectively 
with n,. and y being averaging operators in the corresponding directions 

(ref. 7) (which are defined in the Appends). Then. central difference 
approximations for the derivatives in equation (1) in the current time plane 
at the point ij may be written 
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where hjr and h n are the mesh spacings in the £ and n directions 
respectively. Corresponding formulae apply in the previous time plane except 
that the values of appearing on the right sides of equations (7) are 

evaluated at t* and would therefore be written with asterisks. 


The conventional Crank-Ni col son finite-difference approximation to 
equation (1) at a typical point in the mesh may be written according to 
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where the asterisk denotes known values in the previous time plane and the 
overbar indicates a simple average between the current and the previous time 
plane (as, for example, in equation (6)). It is worthwhile to note that these 
difference approximations are second order accurate in both spatial directions 
and in time. Since the coefficients in equation (8) (Rm, 5m, Tm,...) 
depend in general on the dependent variable u, the difference equations are 
nonlinear and must be solved iteratively; for given estimates of the 
functional coefficients in equation (8), the system of equations is generally 
solved by systematically sweeping the spatial mesh using a Gauss-Seidel or an 
over-relaxation method. In the case when S E 0 in equation (1), the system 
of equation (1) may be written in the form 


u, . = 




fii b ij U i+l,j " C ij U i-l,j d ij U i,j+1 e ij U i,j -If* 


(9) 


at a typical point in the mesh; for equations containing a second order cross 
derivative (and hence having a nonzero S), the right side of equation (9) will 
also involve terms with Ui + i , j + 1 , Ui + i j-l, Uj _i j + ] and , i_i - In a 
point-by-point iterative sweep through the spatial mesh, current estimates of 
the values on the right side of equation (9) are used to define the new 
iterate for U}j at each mesh point. Sweeping of the mesh is continued at 
each time step until convergence is obtained and this typically may require 
the order of 10 to 50 iterations (ref. 1). Generally, convergence occurs 
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rapidly if the system of difference equations is diagonally dominant; this 
implies that the coefficient of the pivotal element must be greater or equal 
in magnitude than the sum of absolute values of the coefficients of the other 
elements in the difference equation. For the system defined by equation (9), 
diagonal dominance occurs if 

l a ijl 2 |bijl + i c ijl + Mijl + l e ijl‘ ^ 

for each point in the mesh. 


Upwind-Downwind Differencing 

It is well known that the central difference approximations (7c) and <7d) 
for the first derivative terms in equation (1) can lead to difference 
equations of the form (9) which are not diagonally dominant. In situations 
where P ^ -i and Qm are large (corresponding to strong local convection 
effects) .the difference equations associated with equations (8) may fail to 
be diagonally dominant over a significant region of the flow field and the 
iterative sweeping of the mesh will fail to converge (refs. 1 and 2). In the 
problem studied by Doligalski and Walker (ref. 1) this phenomenon occurs as a 
viscous boundary- 1 ayer flow near a wall procedes into a strong interaction 
with an outer effectively inviscid flow; in such cases, local strong updrafts 
occur in the boundary layer near the wall. Doligalski and Walker (ref. 1) 
have described an upwind-downwind differencing procedure which is second order 
accurate in both space and time and which always leads to a diagonally 
dominant set of difference equations (see also refs. 8 and 9). The procedure 
only affects the treatment of the terms P3u/3n and Q9n/3£ in equation (1) and 
will be described briefly here to set the stage for the algorithm that will be 
described in the next section. 


Consider first the term P3u/3q. In the conventional Crank-Ni col son 
method, each term in equation (1) is approximated at a point midway between 
the current and previous time planes; a simple average for all terms is then 
carried out along the averaging path labelled (a) in figure 2(a). However, it 
is easily demonstrated that an average along any path through the central point 
yields a second order accurate result. In particular consider the situation 
where Pm > 0 at the typical mesh point and the averaging path labelled (b) 
in figure J 2(a); this path intersects the current time plane at a point midway 
between the points (i,j+l> and (i,j>, and the previous time plane midway 
between the points (i,j) and (i,j-l). Consequently central differences may be 
used for 3U/3n in both the current and previous time planes. For P,j < 0, 
the averaging takes place along the line labelled (c) in figure 2(a). The 
difference approximations for P3u/3n then may be written according to 
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It is easily confirmed that these approximations always act to enhance the 
diagonal dominance of the difference equations; in addition, the 
approximations are second order accurate in both h n and At. A similar 
approach may be adopted for Q 9u/9^ and the averaging paths in this case are 
indicated schematically in figure 2(b); here the average for 9u/9£ takes 
place along either path (d) or path (e) depending on the sign of Oil- The 
difference approximations are J 
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which are second order accurate in h^ and At. 


It Is convenient to write these difference approximations in a more 
compact form by first Introducing the enlargement operators (ref. 7) in the 
£ and n directions defined by 


E y> u(£,n) 
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respectively. Now define operators x?j and Xp according to 
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sgn(PM) = { J 
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It follows that equation (11) may be written 
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In a similar manner, define operators x| and x£ according to 
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(17a) 


xl = 


sgn(Q^/2) 
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and It follows that equation (12) may be written in the form 


x^= 


0 . . — = 

35 

The present algorithm will now be considered. 
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An Upwind-Downwind ADI Method 

Suppose now that the conventional Crank-Nicolson approach is used for all 
terms in equation (1), except P9u/9n and Q9u/9E, where the averaging of the 
spatial partial derivatives takes place along the line labelled (a) in figure 
2(a); the first derivative terms are approximated using equations 
(16) and (18). The resulting finite difference approximation to equation (1) 
may be written 
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(19) 

This equation defines a difference equation at the typical mesh point, which 
can be utilized to carry out a point-by-point sweep through the mesh at each 
time step. Experience with this algorithm (refs. 1 and 2), for the unsteady 
boundary- 1 ayer equations in Eulerian coordinates, indicates that the method 
produces results which are essentially the same as those produced with the 
conventional Crank-Nicolson method for well behaved unsteady boundary-layer 
flows. However, in situations where strong convective effects occur, the 
standard Crank-Nicolson method fails to converge (refs. 1 and 2) as the 
difference equations begin to lose the diagonal dominance property; in such 
cases the upwind-downwind differencing scneme continued to produce converged 
results at each time step and the method could be used well beyond the point 
in time where the conventional Crank-Nicolson method failed (refs. 1 and 2). 
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One disadvantage of the point-by-point iterative methods is that a 
significant number of iterations may be required to obtain convergence at each 
time step. An attractive alernative is an ADI algorithm of the type described 
by Beam and Warming (ref. 5). In this approach, the operators on the left 
side of equation (19) are factored into operators in the £ and n 
directions. In the case of a second order cross derivative, the factorization 
is not easily accomplished and consequently the difference expressions 
resulting from this term are taken to the right side of the equation. A 
factored form of equation (19) may be obtained according to 
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Note that equations (19) and (20) are not completely equivalent and differ by 
terms 0((a At) 2 ); however, the temporal truncation error associated with the 
original Crank-NI col son method is also 0((At) 2 ) and thus equation (20) may be 
regarded as a second order accurate difference approximation for equation (1). 


The set of difference equations given by equation (20) may be solved as a 
sequence of diagonally dominant tridiagonal matrix problems in the following 
manner. Let u^j be an intermediate dependent variable defined by 


U ij " ^ ( R ij 5 C + h ^ij x j) U ij = D 1j’ 


(23) 


it then follows that equation (20) becomes 
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To initiate a solution for Uii in the current time plane, equation (23) is 
first solved along all lines of constant n using a direct method of solution 
for triadiagonal matrix problems (for example, the Thomas algorithm); this 
sweep through the mesh defines current estimates of the intermediate variable 
u n at all internal spatial mesh points Note that boundary conditions for 
Dm along the maximum and minimum values of n for this sweep may be obtained 
using equation (24). In the second phase of the procedure, the computed 
values of u^j are used on the right side of equation (23) which now defines 
a sequence of tridiagonal matrix problems for Ujj along lines of constant £. 
Again this sequence of problems is solved by a direct method for each value 
of £ in the mesh. The net result of this procedure is an estimate for each 
u i -i at each internal mesh point. Note that another version of the algorithm 
isobtained by interchanging the operators in equations (23) and (24) and for 
which a set of tridiagonal problems is first solved along lines of constant 
this would be followed by solving a set of tridiagonal problems along lines 
of constant n- 

It should be noted that it is possible to apply this algorithm in a 
noniterative manner by selecting the time step to be small enough so that 
iteration is not necessary. However, in flow problems in which rapid changes 
begin to develop in the flow field, it is prudent to use a limited amount of 
iteration at each timestep. Iteration is indicated since the coefficients 
Rii, Oii, Tii, and ?ij in equations (23) and (24) are averages between the 
currentand previous time planes and consequently are implicit functions of 
the dependent variable u. Futhermore, the right side of equation (23), which 
is defined in equation (22), contains values of uij explicitly, in terms 
associated with the second order cross derivative; consequently the value of 
Dij at each mesh point should be recomputed and at least one more iteration 
carried out. 


DISCUSSION 

In this study, an upwind-downwind ADI method has been developed. In ADI 
methods, the two-dimensional spatial difference operators are factored into 
two sequential sets of tridiagonal matrix problems, one in each of the 
coordinate directions. These types of methods have proved to be very 
efficient due to the fact that (1) direct solvers for tridiagonal matrix 
problems are relatively fast and accurate, and (2) the process of solving for 
a line of information at once appears to be more effective at communicating 
boundary information to internal mesh points as opposed to systematic 
point-by-point sweeping of the mesh. 

Conventional ADI methods of the type given by Beam and Warming (ref. 5) 
may also be applied to equations like equation (1). Such methods employ 
standard central differences for the spatial derivatives and are generally 
used in a noniterative mode. The present method, as well as the explicit Beam 
and Warming (ref. 5) algorithm, has recently been used to compute the unsteady 
boundary-layer development due to a vortex convected above a wall (ref. 6); at 
a certain stage in this flow, the boundary layer begins to thicken very 
rapidly and strong convective effects develop locally. Indeed this is the 
type of boundary-layer flow which rapidly evolves into an vi scous-i nvi scid 
interaction with the outer inviscid flow in the form of a boundary-layer 
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eruption. In the early stages of the motion, both methods were able to 
compute the flow evolution very efficiently. However, once strong convection 
effects began to develop locally, the Beam and Warming (ref. 7) algorithm was 
unable to successfully track the flow evolution; this failure is believed due 
to a lack of diagonal dominance which eventually leads to failure of the 
method. 

The present method was also compared to the solution process of 
systematic point-by-point iterative sweeping of the mesh. Here again the 
comparison dramatically favored the present method for the example problem 
(ref. 6). Time steps ten times larger could be taken using the present scheme 
and at the same time only about one tenth the iterations were required at each 
time step. Finally, as the convective effects strengthened it proved 
impossible to get converged solution at all, using the point-by-point 
iteration. 

Lastly, it is worthwhile to note that although the present method may be 
used in a noniterative manner, use of iteration at each time step is 
considered prudent and is recommended. In the unsteady boundary- layer problem 
studied in reference 6, a maximum number of iterations per time step was 
assigned a priori ; if the maximum was exceeded for any given time step, the 
calculation was restarted from the previous time plane with a smaller time 
step. In this manner, the number of iterations per time step gives an 
indication of the need to reduce (or increase) the time step. Ideally the 
time step should be such that on the order of two or three iterations only are 
required each time step. If very little change occurs in the solution over 
each iteration, a larger time step can be used. 
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APPENDIX 


In this Appendix, the central difference operator notation used in this 
study is explicitly defined for completeness. Here (£,n> are independent 
spatial coordinates and the uniform mesh spacing in each direction Is noted by 

h and h respectively. The enlargement operators In the £ and n 
t n 

directions are defined by 

E^uCE.n) = u(E - ph^,n>, <A.1) 

E q u(E,n> = u<E,n + qh >, <A.2) 

n n 


where p and q are constants; these operators simply shift the appropriate 
argument of the function. The central difference operators and 6^ are 

defined by 
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5 = 


i! 1 - e 1/3 



(A. 3) 


and consequently 


S u^(^,n) = u(E 

S u(E,n) = u(£,n + 
n 


+ n^/2,n> - u<£ -^h /2,n>, 
h / 2) - u<E,n - h /2) , 


(A. 4) 
(A. 5) 


corresponding to a central difference with one spatial variable held fixed. 
Finally the averaging operators and p n are defined by 


Pc - 


1 c 1/j p-* 

2 e e + h 


P„ = 


± E + E 
2 n n 


(A. 6) 


(A. 7) 


and therefore 


p^.u(E,n> = ^ju(E + - h^/2,n) 

p^uCE.n) = | ju(E .n + h^/2) + u(E,n - b^/2) 


(A. 8) 
(A. 9) 


corresponding to an average with one spatial variable held fixed. 
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FIGURE 1. - GRID STRUCTURE AND NOTATION FOR CONVENTIONAL CRANK-NICOLSON METHOD. 




(b) G du/d£. 

FIGURE 2. - AVERAGING PATH FOR CONVECTIVE DERIVATIVE TERMS. 
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